using ControlSystemsBase, LinearAlgebra, Plots, Statistics, LaTeXStrings
using Plots.PlotMeasures

First, we start by a simple LQR control design of the following linear system:

Δ = 0.1
A = [0.98 Δ; 0 0.95];
B = [0; Δ;;];
C = I;
D = 0;
sys = ss(A,B,C,D,Δ)

V = [0.5;;]
W = diagm([0.2, 0.05])

Q = I
R = 2I

K_LQR = -lqr(sys, Q, R)

display(K_LQR)

1×2 Matrix{Float64}:
 -0.441338  -0.788109

We can also achieve the same optimal controller using the controllability-type Gramian, i.e., via LMIs in an SDP (Problem~2 in the paper):

using JuMP
using MosekTools, SCS

model = Model(Mosek.Optimizer);
rₓ = size(A)[1]
rᵤ = size(B)[2]
@variable(model, Σ[i=1:rₓ, j=1:rₓ], PSD)
@variable(model, Z₀[i=1:rᵤ, j=1:rᵤ], PSD)
@variable(model, L[i=1:rᵤ, j=1:rₓ])
@objective(model, Min, tr(Q*Σ) + tr(R*Z₀))
@constraint(model, Σ-1e-5*I >= 0, PSDCone())
@constraint(model, [Z₀ L; L' Σ] >= 0, PSDCone())
# Lyupanov
@constraint(model, Σ - (A*Σ*A' + B*L*A' +
 A*L'*B' +  B*Z₀*B' + W + B*V*B') >= 0, PSDCone());
@constraint(model, Σ - (A*Σ*A' + B*L*A' +
 A*L'*B' +  B*Z₀*B' + W + B*V*B') <= 0, PSDCone());

set_silent(model)
optimize!(model);

L, Σ, Z₀ = value.(L), value.(Σ), value.(Z₀)
K_LMI = L * inv(Σ)

1×2 Matrix{Float64}:
 -0.44132  -0.78812

model = Model(Mosek.Optimizer);
rₓ = size(A)[1]
rᵤ = size(B)[2]
@variable(model, Σ[i=1:rₓ, j=1:rₓ], PSD)
@variable(model, Z₀[i=1:rᵤ, j=1:rᵤ], PSD)
@variable(model, L[i=1:rᵤ, j=1:rₓ])
@objective(model, Min, tr(Q*Σ) + tr(R*Z₀))
@constraint(model, Σ-1e-5*I >= 0, PSDCone())
@constraint(model, [Z₀ L; L' Σ] >= 0, PSDCone())
# Lyupanov
@constraint(model, [(Σ-W-B*V*B') A*Σ+B*L; (A*Σ+B*L)' Σ]>=0, PSDCone());

set_silent(model)
optimize!(model);

L, Σ, Z₀ = value.(L), value.(Σ), value.(Z₀)
K_LMI = L * inv(Σ)

1×2 Matrix{Float64}:
 -0.441352  -0.788199

The error between the two control gains:

norm(K_LMI - K_LQR, 2) # the two norm of the difference.

2.032880629709027e-5

We now seek to do the previous step again, but with $(A,B)$ being identified from data.

First, we use the double integrator model above, to generate noisy data, under a stabilizing controller $K_0$ and a persistently exciting additive input $v_k$. That is, $u_k = K_0 x_k + v_k$.

function f(x::Vector{Float64}, u::Vector{Float64})
    x⁺ = zero.(x);
    x⁺ = A * x + B * u + sqrt(W) * randn(size(x⁺))
    return x⁺
end

f (generic function with 1 method)

function f_dataGen_CL(x₀::Vector{Float64}, K::Matrix{Float64}, N::Int)
    X = zeros(Float64, size(K)[2], N+1)
    U = zeros(Float64, size(K)[1], N)
    X[:,1] = x₀
    for j=1:N
        U[:,j] = K * X[:,j] + sqrt(V) * randn(rᵤ,)
        X[:,j+1] = f(X[:,j], U[:,j])
    end
    return X[:,1:end-1], U
end

f_dataGen_CL (generic function with 1 method)

N = 2000
x₀ = zeros(Float64, 2,)
K₀ = [-0.2  -9.0]

X, U = f_dataGen_CL(x₀, K₀, N)
display(abs.(eigen(A + B * K₀).values))

plt = plot(
    framestyle = :box,
    yguidefontsize = 8,
    xguidefontsize = 8,
    xtickfontsize = 8,
    ytickfontsize = 8,
    ylabel = L"$x^2_k$",
    xlabel = L"$x^1_k$",
    palette = :seaborn_muted,
    legend=:none,
    foreground_color_legend = nothing,
    legendfontsize=5,
    fontfamily = "Computer Modern",
    markersize = 2,
)
plot_data_exp = scatter!(plt, X[1,:], X[2,:], markersize = 3) 
#[:none, :auto, :circle, :rect, :star5, :diamond, :hexagon, :cross, :xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon, :heptagon, :octagon, :star4, :star6, :star7, :star8, :vline, :hline, :+, :x].

2-element Vector{Float64}:
 0.05215553368329004
 0.97784446631671

We use the above input and state data, $U$ and $X$, respectively, to fit a linear model and design a new controller based on this identified model, in a certainty equivalence fashion.

AB = X[:,2:end] / [X[:,1:end-1]; U[:,1:end-1]]
 = AB[:,1:rₓ]
B̂ = AB[:,rₓ+1:end]
sys_estimate = ss(Â,B̂,I,0,Δ)

K_CE = -lqr(sys_estimate,Q,R)

1×2 Matrix{Float64}:
 -0.515261  -0.712289

We now use this certainty equivalence control law to control the true system.

X_CE, _ = f_dataGen_CL(x₀, K_CE, N)
plot_CE = scatter(plot_data_exp, X_CE[1,:], X_CE[2,:], markersize = 3, markershape = :circle)

Notice that the two distributions do not share the same support exactly. That is, the new system visits states that were not experienced during the experiment and learning phase. These regions of the state-space may contain nonlinearities, discontinuities, unpredicted safety violations. The responsibility of data-conforming control is to prevent such sudden generalization beyond learning data, and to enforce that the new system's distribution is close (under some norm) to the data distribution.

function solveProb5(γ_prime)
    Σ_data = 1/(N) * X * X'
    model = Model(Mosek.Optimizer);
    rₓ = size(A)[1]
    rᵤ = size(B)[2]

    @variable(model, Σ[i=1:rₓ, j=1:rₓ], PSD)
    @variable(model, Z_F[i=1:rₓ, j=1:rₓ], PSD)
    @variable(model, Z₀[i=1:rᵤ, j=1:rᵤ], PSD)
    @variable(model, L[i=1:rᵤ, j=1:rₓ])

    @objective(model, Min, tr(Q*Σ) + tr(R*Z₀) + γ_prime * tr(Z_F))

    @constraint(model, Σ - (1e-6I) >= 0, PSDCone())
    @constraint(model, [Z₀ L; L' Σ] >= 0, PSDCone())
    # Lyupanov
    #@constraint(model, Σ - (Â*Σ*Â' + B̂*L*Â' + Â*L'*B̂' +  B̂*Z₀*B̂' + W + B̂*V*B̂') >= 0, PSDCone());
    #@constraint(model, Σ - (Â*Σ*Â' + B̂*L*Â' + Â*L'*B̂' +  B̂*Z₀*B̂' + W + B̂*V*B̂') <= 0, PSDCone());
    @constraint(model, [(Σ-W-B̂*V*B̂') Â*Σ+B̂*L; (Â*Σ+B̂*L)' Σ]>=0, PSDCone());
    @constraint(model, [Z_F (Σ-Σ_data); (Σ-Σ_data) I] >= 0, PSDCone());
    #set_attribute(model, "INTPNT_CO_TOL_DFEAS", 1e-10)
    set_silent(model)
    optimize!(model);

    L, Σ, Z₀ = value.(L), value.(Σ), value.(Z₀)
    K_prob5 = L * inv(Σ)
    
    return K_prob5, Σ
end

solveProb5 (generic function with 1 method)

We solve Problem~5 with $\gamma$ set to $0$, expecting the same behavior as the solution of the CE-LQR problem.

γ_prime = 5
K_prob5, _ = solveProb5(γ_prime)
display(K_prob5)
X_prob5_0, _ = f_dataGen_CL(x₀, K_prob5, N)
plot_prob5_0 = scatter(plot_data_exp, X_prob5_0[1,:], X_prob5_0[2,:], markersize = 2)

1×2 Matrix{Float64}:
 -0.0797656  -0.5735

We now start to increase $\gamma$ and notice the difference.

γ_prime = 20
K_prob5, _ = solveProb5(γ_prime)
X_prob5_1, _ = f_dataGen_CL(x₀, K_prob5, N)
plot_prob5_1 = scatter(plot_data_exp, X_prob5_1[1,:], X_prob5_1[2,:], markersize = 2)

γ_prime = 100
K_prob5, _ = solveProb5(γ_prime)
X_prob5_2, _ = f_dataGen_CL(x₀, K_prob5, N)
plot_prob5_20 = scatter(plot_data_exp, X_prob5_2[1,:], X_prob5_2[2,:], markersize = 2, order = 1)

Plotting all of the above figures in the same figure and with the same axes scaling:

x_minimum = minimum([minimum(X[1,:]), minimum(X_prob5_0[1,:]), minimum(X_prob5_1[1,:]), minimum(X_prob5_2[1,:])])
x_maximum = maximum([maximum(X[1,:]), maximum(X_prob5_0[1,:]), maximum(X_prob5_1[1,:]), maximum(X_prob5_2[1,:])])
x_limits = (x_minimum, x_maximum)

y_minimum = minimum([minimum(X[2,:]), minimum(X_prob5_0[2,:]), minimum(X_prob5_1[2,:]), minimum(X_prob5_2[2,:])])
y_maximum = maximum([maximum(X[2,:]), maximum(X_prob5_0[2,:]), maximum(X_prob5_1[2,:]), maximum(X_prob5_2[2,:])])
y_limits = (y_minimum, y_maximum)

plt = plot(
    framestyle = :box,
    yguidefontsize = 6,
    xguidefontsize = 6,
    xtickfontsize = 6,
    ytickfontsize = 6,
    ylabel = L"$x^2_k$",
    xlabel = L"$x^1_k$",
    palette = :seaborn_muted,
    legend=:none,
    foreground_color_legend = nothing,
    legendfontsize=5,
    fontfamily = "Computer Modern",
    markersize = 2,
    titlefontsize = 6,
)
MarkerSize = 2
plot_CE = scatter(plt, X_CE[1,:], X_CE[2,:], markersize = MarkerSize, markershape = :xcross, bottom_margin = [-2mm 0mm], title = "(a)")
plot_prob5_0 = scatter(plt, X_prob5_0[1,:], X_prob5_0[2,:], markersize = MarkerSize, markershape = :xcross, bottom_margin = [-2mm 0mm], title = "(b)")
plot_prob5_1 = scatter(plt, X_prob5_1[1,:], X_prob5_1[2,:], markersize = MarkerSize, markershape = :xcross, bottom_margin = [-2mm 0mm], title = "(c)")
plot_prob5_20 = scatter(plt, X_prob5_2[1,:], X_prob5_2[2,:], markersize = MarkerSize, markershape = :xcross, bottom_margin = [-2mm 0mm], title = "(d)")
plots_list = [plot_CE, plot_prob5_0, plot_prob5_1, plot_prob5_20]
for i = 1:length(plots_list)
    scatter!(plots_list[i], X[1,:], X[2,:], xlims = x_limits, ylims = y_limits, xlabel = "", ylabel = "", 
    xticks = [ceil(x_minimum), 0, floor(x_maximum)], yticks = [ceil(y_minimum), 0, floor(y_maximum)], 
    markersize = MarkerSize, markerstrokewidth=1, markershape = :rect)
end

width = 320
height = 360
Fig1 = plot(plot_CE, plot_prob5_0, plot_prob5_1, plot_prob5_20, layout = (2,2), size = (width, height))

Fig1
savefig("figs/Prob5.pdf")

"/Users/mohammad/Documents/GitHub/Data-Conforming-Control/DataConformingControl/figs/Prob5.pdf"

Now, we try to solve the corresponding Problem~6 of this problem above.

function solveProb6(γ)
    rₓ = size(A)[1]
    rᵤ = size(B)[2]
    D = [X;U]
    Σ_data = 1/(N) * X * X'
    Γ_data = 1/(N) * D * D'
    H_data = Γ_data[1:rₓ,rₓ+1:end]
    Υ = cat(V, zeros(rₓ,rₓ); dims=(1,2))
    model = Model(Mosek.Optimizer);
    

    @variable(model, Σ[i=1:rₓ, j=1:rₓ], PSD)
    @variable(model, Z₀[i=1:rᵤ, j=1:rᵤ], PSD)
    @variable(model, L[i=1:rᵤ, j=1:rₓ])
    @variable(model, Z₁[i=1:rₓ+rᵤ, j=1:rₓ+rᵤ])
    @variable(model, Z₂[i=1:rᵤ, j=1:rᵤ])
    @variable(model, Z₃[i=1:rₓ, j=1:rₓ])

    @objective(model, Min, tr(Q*Σ) + tr(R*Z₀) 
    + γ * (tr(inv(Γ_data) * Z₁) + tr(inv(V) * Z₂) + tr(Σ_data * Z₃)))

    @constraint(model, Σ - (1e-6I) >= 0, PSDCone())
    @constraint(model, [Z₀ L; L' Σ] >= 0, PSDCone())
    # Lyupanov
    #@constraint(model, Σ - (Â*Σ*Â' + B̂*L*Â' + Â*L'*B̂' +  B̂*Z₀*B̂' + W + B̂*V*B̂') >= 0, PSDCone());
    #@constraint(model, Σ - (Â*Σ*Â' + B̂*L*Â' + Â*L'*B̂' +  B̂*Z₀*B̂' + W + B̂*V*B̂') <= 0, PSDCone());
    @constraint(model, [(Σ-W-B̂*V*B̂') Â*Σ+B̂*L; (Â*Σ+B̂*L)' Σ]>=0, PSDCone());
    # The Z's
    @constraint(model, [Z₃ I; I Σ] >= 0, PSDCone())
    @constraint(model, [Z₂ (L-H_data'*inv(Σ_data)*Σ); (L-H_data'*inv(Σ_data)*Σ)' Σ] >= 0, PSDCone())
    @constraint(model, [(Z₁-Υ) cat(Σ_data, L; dims=1); cat(Σ_data, L; dims=1)' Σ] >= 0, PSDCone())
    #set_attribute(model, "INTPNT_CO_TOL_DFEAS", 1e-10)
    set_silent(model)
    optimize!(model);

    L, Σ, Z₀ = value.(L), value.(Σ), value.(Z₀)
    K_prob6 = L * inv(Σ)
    
    return K_prob6, Σ
end

solveProb6 (generic function with 1 method)

γ = 100
K_prob6, _ = solveProb6(γ)
display(K_prob6)
X_prob6, _ = f_dataGen_CL(x₀, K_prob6, N)
plot_prob6 = scatter(plot_data_exp, X_prob6[1,:], X_prob6[2,:], markersize = 2, order = 1)

1×2 Matrix{Float64}:
 -0.179531  -8.30274

Now, to make things more interesting, we add a nonlinearity to the linear system (the function $f$) that will render a non-data-conforming control unstable.

function f(x::Vector{Float64}, u::Vector{Float64})
    x⁺ = zero.(x);
    x⁺ = A * x + B * u + sqrt(W) * randn(size(x⁺))
    x⁺[1] += 1/9 * x[2] ^ 2
    return x⁺
end

f (generic function with 1 method)

We now generate data using the simulation experiment with the initial control gain $K_0$.

N = 2000
x₀ = zeros(Float64, 2,)
K₀ = [-0.2  -9.0]

X, U = f_dataGen_CL(x₀, K₀, N)
x₀ = X[:,end];

We fit a new model using the new data, and solve Problem~3.

AB = X[:,2:end] / [X[:,1:end-1]; U[:,1:end-1]]
 = AB[:,1:rₓ]
B̂ = AB[:,rₓ+1:end]
sys_estimate = ss(Â,B̂,I,0,Δ)

K_CE = -lqr(sys_estimate,Q,R)
X_CE, _ = f_dataGen_CL(x₀, K_CE, N);

We also solve Problem~6 with multiple values of $\gamma$:

γ_prime = 5.0
K_prob5, _ = solveProb5(γ_prime)
X_prob5_0, _ = f_dataGen_CL(x₀, K_prob5, N);

γ_prime = 20
K_prob5, _ = solveProb5(γ_prime)
X_prob5_1, _ = f_dataGen_CL(x₀, K_prob5, N);

γ_prime = 100
K_prob5, _ = solveProb5(γ_prime)
X_prob5_2, _ = f_dataGen_CL(x₀, K_prob5, N);

x_minimum = minimum(X[1,:]) - 8
x_maximum = maximum(X[1,:]) + 8

plt = plot(
    framestyle = :box,
    yguidefontsize = 6,
    xguidefontsize = 6,
    xtickfontsize = 6,
    ytickfontsize = 6,
    ylims = (x_minimum, x_maximum),
    palette = :seaborn_muted,
    legend=:none,
    foreground_color_legend = nothing,
    legendfontsize=5,
    fontfamily = "Computer Modern",
    markersize = 2,
    titlefontsize = 6,
)
MarkerSize = 2
plot_CE = plot(plt, [X[1,:];X_CE[1,:]], bottom_margin = [-2mm 0mm], ylims = (x_minimum, x_maximum), xticks = :none, ylabel = "(1)")
plot_prob5_0 = plot(plt, [X[1,:];X_prob5_0[1,:]], bottom_margin = [-2mm 0mm], ylims = (x_minimum, x_maximum), xticks = :none, ylabel = "(2)")
plot_prob5_1 = plot(plt, [X[1,:];X_prob5_1[1,:]], bottom_margin = [-2mm 0mm], ylims = (x_minimum, x_maximum), xticks = :none, ylabel = "(3)")
plot_prob5_20 = plot(plt, [X[1,:];X_prob5_2[1,:]], bottom_margin = [-2mm 0mm], ylims = (x_minimum, x_maximum), xlabel = L"$k$", ylabel = "(4)")
plots_list = [plot_CE, plot_prob5_0, plot_prob5_1, plot_prob5_20]
for i = 1:length(plots_list)
    #plot!(plots_list[i], X[1,:], ylims = (x_minimum, x_maximum))
end

width = 320
height = 320
Fig2 = plot(plot_CE, plot_prob5_0, plot_prob5_1, plot_prob5_20, layout = (4,1), size = (width, height))

Fig2
savefig("figs/Prob5_nonlinear.pdf")

"/Users/mohammad/Documents/GitHub/Data-Conforming-Control/DataConformingControl/figs/Prob5_nonlinear.pdf"

We now repeat the whole procedure above bunch of times to get some statistics about the stability of the above control design procedures.

stability_count = zeros(4,)
num_of_simulations = 1000
for j = 1:num_of_simulations
stability_condition = false
while ~stability_condition
    x₀ = zeros(Float64, 2,)
    K₀ = [-0.2  -9.0]
    X, U = f_dataGen_CL(x₀, K₀, N)
    # check stability (boundedness) of the experimental data:
    maximum_state = maximum(abs.(X))
    if maximum_state >= 50
        stability_condition = false
    else
        stability_condition = true
    end
end
AB = X[:,2:end] / [X[:,1:end-1]; U[:,1:end-1]]
 = AB[:,1:rₓ]
B̂ = AB[:,rₓ+1:end]
sys_estimate = ss(Â,B̂,I,0,Δ)

# CE control
K_CE = -lqr(sys_estimate,Q,R)
X_CE, _ = f_dataGen_CL(x₀, K_CE, N)
maximum_state = maximum(abs.(X_CE))
    if maximum_state <= 50
        stability_count[1] += 1
    end
# Data conforming
γ_prime = 5.0
K_prob5, _ = solveProb5(γ_prime)
X_prob5_0, _ = f_dataGen_CL(x₀, K_prob5, N);
maximum_state = maximum(abs.(X_prob5_0))
    if maximum_state <= 50
        stability_count[2] += 1
    end

γ_prime = 20
K_prob5, _ = solveProb5(γ_prime)
X_prob5_1, _ = f_dataGen_CL(x₀, K_prob5, N);
maximum_state = maximum(abs.(X_prob5_1))
    if maximum_state <= 50
        stability_count[3] += 1
    end
γ_prime = 100
K_prob5, _ = solveProb5(γ_prime)
X_prob5_2, _ = f_dataGen_CL(x₀, K_prob5, N);
maximum_state = maximum(abs.(X_prob5_2))
    if maximum_state <= 50
        stability_count[4] += 1
    end

end

println("The following are the percentages of stable simulations by each controller:")
display(stability_count./num_of_simulations)

4-element Vector{Float64}:
 0.26
 0.996
 0.999
 1.0

The following are the percentages of stable simulations by each controller:

We fit a new model using the new data, and solve Problem~3.

We also solve Problem~6 with multiple values of $\gamma$:

State-input joint data-conforming:

Suppose now that the state dynamics contain cross state-input terms, for instance:

function f(x::Vector{Float64}, u::Vector{Float64})
    x⁺ = zero.(x);
    x⁺ = A * x + B * u + sqrt(W) * randn(size(x⁺))
    x⁺[1] += 1/9 * x[2] ^ 2
    x⁺[2] +=  1/9*tanh(1*x[1]) * u[1]
    return x⁺
end

f (generic function with 1 method)

We now generate data using the simulation experiment with the initial control gain $K_0$.

N = 2000
x₀ = zeros(Float64, 2,)
K₀ = [-0.2  -9.0]

X, U = f_dataGen_CL(x₀, K₀, N)
x₀ = X[:,end];
scatter(X[1,:], X[2,:])

We fit a new model using the new data, and solve Problem~3.

AB = X[:,2:end] / [X[:,1:end-1]; U[:,1:end-1]]
 = AB[:,1:rₓ]
B̂ = AB[:,rₓ+1:end]
sys_estimate = ss(Â,B̂,I,0,Δ)

K_CE = -lqr(sys_estimate,Q,R)
X_CE, _ = f_dataGen_CL(x₀, K_CE, N);

We also solve Problem~6 with multiple values of $\gamma$:

γ_prime = 0.0
K_prob5, _ = solveProb5(γ_prime)
X_prob5_0, _ = f_dataGen_CL(x₀, K_prob5, N);

γ_prime = 20
K_prob5, _ = solveProb5(γ_prime)
X_prob5_1, _ = f_dataGen_CL(x₀, K_prob5, N);

γ_prime = 100
K_prob5, _ = solveProb5(γ_prime)
X_prob5_2, _ = f_dataGen_CL(x₀, K_prob5, N);

γ = 10
K_prob6, _ = solveProb6(γ)
X_prob6, _ = f_dataGen_CL(x₀, K_prob6, N)
plot_prob6 = scatter(plot_data_exp, X_prob6[1,:], X_prob6[2,:], markersize = 2, order = 1)

x_minimum = minimum(X[1,:]) - 12
x_maximum = maximum(X[1,:]) + 12

plt = plot(
    framestyle = :box,
    yguidefontsize = 6,
    xguidefontsize = 6,
    xtickfontsize = 6,
    ytickfontsize = 6,
    ylims = (x_minimum, x_maximum),
    palette = :seaborn_muted,
    legend=:bottomleft,
    foreground_color_legend = nothing,
    legendfontsize=6,
    fontfamily = "Computer Modern",
    markersize = 2,
    titlefontsize = 6,
)
MarkerSize = 2
plot_CE = plot(plt, [X[1,:];X_CE[1,:]], bottom_margin = [-2mm 0mm], ylims = (x_minimum, x_maximum), label = "Prob. 3")
plot_prob5_1 = plot(plot_CE, [X[1,:];X_prob5_1[1,:]], bottom_margin = [-2mm 0mm], ylims = (x_minimum, x_maximum), label = L"Prob. 5 ($\gamma'=10$)")
plot_prob5_20 = plot(plot_prob5_1, [X[1,:];X_prob5_2[1,:]], bottom_margin = [-2mm 0mm], ylims = (x_minimum, x_maximum), label = L"Prob. 5 ($\gamma'=100$)")
plot_prob6 = plot(plot_prob5_20, [X[1,:];X_prob6[1,:]], bottom_margin = [0mm 0mm], 
ylims = (x_minimum, x_maximum), xlabel = L"$k$", label = L"Prob. 6 ($\gamma=10$)", linestyle = :dashdot)


width = 320
height = 160
Fig3 = plot(plot_prob6, size = (width, height))

Fig3
savefig("figs/Prob5_bilinear.pdf")

"/Users/mohammad/Documents/GitHub/Data-Conforming-Control/DataConformingControl/figs/Prob5_bilinear.pdf"

Similar to the previous example, we run a bunch of simulations and count the number of stable simulations offered by each control:

stability_count = zeros(4,)
num_of_simulations = 1000
for j = 1:num_of_simulations
stability_condition = false
while ~stability_condition
    x₀ = zeros(Float64, 2,)
    K₀ = [-0.2  -9.0]
    X, U = f_dataGen_CL(x₀, K₀, N)
    # check stability (boundedness) of the experimental data:
    maximum_state = maximum(abs.(X))
    if maximum_state >= 50
        stability_condition = false
    else
        stability_condition = true
    end
end
AB = X[:,2:end] / [X[:,1:end-1]; U[:,1:end-1]]
 = AB[:,1:rₓ]
B̂ = AB[:,rₓ+1:end]
sys_estimate = ss(Â,B̂,I,0,Δ)

# CE
γ_prime = 0.0
K_prob5, _ = solveProb5(γ_prime)
X_prob5_0, _ = f_dataGen_CL(x₀, K_prob5, N);
maximum_state = maximum(abs.(X_prob5_0))
    if maximum_state <= 50
        stability_count[1] += 1
    end

# State-data conforming
γ_prime = 20
K_prob5, _ = solveProb5(γ_prime)
X_prob5_1, _ = f_dataGen_CL(x₀, K_prob5, N);
maximum_state = maximum(abs.(X_prob5_1))
    if maximum_state <= 50
        stability_count[2] += 1
    end

γ_prime = 100
K_prob5, _ = solveProb5(γ_prime)
X_prob5_2, _ = f_dataGen_CL(x₀, K_prob5, N);
maximum_state = maximum(abs.(X_prob5_2))
    if maximum_state <= 50
        stability_count[3] += 1
    end

# Joint state-input data-conforming
γ = 10
K_prob6, _ = solveProb6(γ)
X_prob6, _ = f_dataGen_CL(x₀, K_prob6, N)
maximum_state = maximum(abs.(X_prob6))
    if maximum_state <= 50
        stability_count[4] += 1
    end

end

println("The following are the percentages of stable simulations by each controller:")
display(stability_count./num_of_simulations)

4-element Vector{Float64}:
 0.078
 0.264
 0.256
 0.981

The following are the percentages of stable simulations by each controller: